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Gravimetric Tidal Loading Computed from Integrated 

Green’s Functions 

Clyde C. Goad 
National Geodetic Survey 
National Ocean Survey, NOAA 
Rockville, Md. 20852 


ABSTRACT, The usual method of predicting the 
effects of ocean tides on geodetic measure- 
ments is to use impulse response functions 
(called Green’s functions) by convolving them 
with the desired ocean tide model. Because 
ocean tide representations are usually ex- 
pressed as areas or cells of constant ampli- 
tude and phase, it has been found that the 
integrals of Green's functions are more 
desirable for use with tidal loading 
calculations . 

Predictions are presented of the loading 
effects on Earth-tide gravimeter measurements 
using a global M2 ocean tide model developed 
by E. W. Schwiderski. Souriau's calculations 
for loading effects on tidal gravity data 
available for western Europe are confirmed to 
an accuracy of ±0.2 microgal for most sites. 
Other results are available from California, 
Australia, and Japan. 


INTRODUCTION 


For several decades prediction of the response of the Earth to 
variable loads has captured the interest of geophysicists, ocean- 
ographers, and geodesists. During the past decade, improvements 
in Earth modeling enabled researchers to determine realistic 
response functions given the elastic properties of the Earth. 
It is only natural to use these tools to predict the response of 
the Earth to the loading of the ocean tides. In this paper an 
alternate form has been developed for the response functions . 
This method requires that the ocean tide height be given in 
terms of areas of constant amplitude and phase, the way global 
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tide models are usually computed. When such representations of 
the tide height are available over any cell or area, the ampli- 
tude times the sine or cosine of phase is constant and can be 
taken outside the integral. Then the integral of the impulse 
response function (or Green's function) remains to be evaluated. 
This integral is stable (by removal of one singularity) and en- 
ables one easily to use any set of load deformation coefficients 
in the study of Earth gravity or displacement response. Previ- 
ously, one had to use published Green's functions (Farrell 1972) 
or compute the functions, requiring the evaluation of infinite 
series which do not converge when the angular argument is small. 
The technique presented here also allows one properly to account 
for the distance of the instrument above sea level. Of course, 
this method is not restricted to global representations. It 
also can be used in any region where the tide is represented by 
areas of constant amplitude and phase. 

The results presented in this paper are based on the global 
representation of the M2 tide developed by E. W* Schwiderski of 
the Naval Surface Weapons Center, Dahlgren, Va. This model is 
constrained to agree with coastal and island tide data, and is 
used here to predict the effect of ocean loading on gravity data 
taken at several areas of the Earth, The results look promising. 
For example, in Australia, the tidal perturbations caused by the 
oceans appear to be predicted to the 0.5-microgal level, 
slightly better than the 1-microgal level reported by Breteger 
and Mather (1978). When the tidal gravimeter stations are more 
than 1° from the land-water boundaries, excellent agreement is 
demonstrated with Souriau (1979) who used 0.25° grids for nearby 
seas when correcting many tidal gravity observations taken in 
western Europe. 
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RESPONSE FUNCTIONS 
Loading Potential 

For the special case of the ocean tides, let the mass distri- 
bution be represented by a constant density layer with varying 
height, h, covering a large sphere of radius, a. The gravita- 
tional potential becomes 


U* 


Gpa 2 



h(6 , a ) sin 6 do da 

■ ■ • 

Va s :• r 2 ‘ - 2ar cos 0 


( 1 ) 


The quantities 0 and a are the central angle and azimuth, r iz 
the distance from t;he center of the Earth, G is the gravitation- 
al constant, and p is the density of sea water. Noting the 
presence of the generating function for the Legendre polynomials 
P n (cos 0), the integral becomes 


= Gpa 2 E jj h(0,a) P n (cos 0) sin 0 d0 da (2a) 


for the solution outside the sphere (exterior), and 


= Gpa 2 E I J h(0,a) P n (cos 0) ■ r , sin 0 d0 da (2b) 

for the interior solution. The solution at the spherically 
coated surface is obtained from eq. (1) by letting r = a. How- 
ever, measurements with Earth tidal gravimeters normally take 
place above sea level, and thus eq. (1) or (2a) must be used. 
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Newtonian Tidal Attraction ^valuation 

Several numerical solutions of the Laplace tidal equations 
have recently become available. Previously, solutions were 
given by cotidal charts which showed contours of constant ampli- 
tude and phase, or low degree spherical harmonic expansions. 
Now the solutions are almost always provided in gridded form 
where the tidal amplitude and phase are given as constants over 
small areas (e.g., 1° geographic squares). Let A^ and be the 
amplitude and Greenwich phase, respectively, over the i-th 
region. Then the tide height, h^(t), at time t is given by 

h i<t) = cos (n • p + 0 .) (3) 

where ii is the coefficient vector for a given constituent and p 
is a vector of six astronomical angles (Cartwright and Tayler 
1971). Decomposing eq. (3) with multiple angle identities 
yields the sinusoidal and cosinusoidal terms 

h^(t) = h? cos (n • p) - h? sin (n • p). (4) 

Substituting eq. (4) into eq. (1), and assuming that the h? and 
h® are constant over the i-th sector bounded by azimuthal angles 
and 0 f 2 ^ and central angles and 0 2l . , yields a rather 

simple representation of the direct loading potential 

U* (a 2 i " a ii) £h?cos(n • p) - h? sin(n • p) 

(5) 

• £,/a 2 + r 2 - 2ar cos 0 

Differentiation of eq. (5), or either 2a or 2b, with respect to 
r yields gravity above or below the tidal layer. The 
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mathematical model of an infinitesimally thin layer covering the 
surface of a sphere exhibits a discontinuity in gravity as one 
crosses the boundary of the tidal sheet. Thus gravity readings 
made very close to coasts or on islands can be significantly per- 
turbed. (See Pekeris (1978) for a further discussion of the 
problem. ) 


Load Deformation Coefficients hi, ki, and 11 

n' n n 

Because the Earth is not a perfectly rigid body, a deformation 
will occur due to the application of the tidal load. The re- 
sponse of the Earth to the load is described by the load defor- 
mation coefficients h^, k^, 1^ (Munk and MacDonald 1960). Let 
represent the n-th degree contribution to the loading poten- 
tial in eq. (2). Then the change in the distance from the 
center of the Earth (or height) is given by h^U^/g. The change 
in the n-th degree potential caused by the redistribution of 
mass is represented by the coefficients k^. The actual poten- 
tial after deformation is given by (l+k^U^. Similarly to h^, 
the 1^ represents horizontal displacements of degree n. 
Numerical values of load deformation coefficients used in this 
study have been taken from Farrell (1972) and Zschau (1978). 

Gravity Change Caused by Deformation 

Since the direct or Newtonian contribution to gravity has been 
given, one must now determine the contribution to gravity result- 
ing from deformation. For the external problem differentiation 
of Sk^n with respect to r gives 


1 

r 


EK (n+l ) O' 
n=o 
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The change in gravity resulting from height changes is given by 

00 

H Eh' U'/g. 
r n n' 3 

n=o 

The combination of these two terms along with the Newtonian con- 
tribution gives the total change in gravity caused by the tidal 
mass layer 


Stifle = ? [(n+l)(l+k A )-2hn U A (6) 

n=o L J 

where the sign is reversed so a downward attraction is positive 
(as is measured by a gravimeter). The term for the attraction 
inside the brace should be used rather than the corresponding 
term for attraction under the tidal sheet given by Favrell 
(1972) and Longman (1963). 

Summation of Series 


As previously discussed, the Newtonian contribution is best 
evaluated by using eq. (5) which avoids the evaluation of an in 
finite sum. The important difference between this approach and 
that given by other investigators is that the integrals of 
Green's functions are computed rather than Green's functions 
themselves. Again noting that the amplitude and phase of the 
ocean tide are constant over limited areas (as was done for the 
Newtonian attraction), the tide height term can be taken outside 
the integral to simplify the process. This technique then re- 
duces to an evaluation of the expression JP n (cos 0) sin 0 d0 
which can be obtained by using recursive expressions. Let 
T n (0)= JP n (cos 0) sin 0 d0 . Then the T's are given by 


v e > 


sin 0 
n(n+l ) 



(cos 


0 ) 


(7) 


7 


where P n ^(cos 6) is the associated Legendre function of degree n 
and order one. The recursive expressions for P n i( cos e ) are 
used to obtain 


T n' e > = HTT cos 0 T n-l< e > * H7T T n-2< 0 >' < 8 > 

The functions T are very desirable in that the infinite sum 

» . . n « 
n-o T n ex ^ 1 ^ lts no singularities as does P n when the central 

angles are small. These are essentially the disk factors that 
Farrell (1972, 1973) used to improve the convergence character- 
istics of Green's functions. Although not required, Hummer's 
method (Farrell 1972) can also be used to facilitate the evalu- 
ation of the infinite series because h^ and nk^ approach con- 
stants as n gets large. The terms involving height aoove 
sea-level, (a/r) n , remain in the infinite sums. These can be 
important especially if the gravimeter is placed in rather high 
locations. Pekeris (1978) has shown that as one approaches the 
surface from above or below the tidal sheet, the infinite sum 
becomes a composition of two terms. One term represents the 
solution in the center of the surface or boundary, and the other 
is a delta function accounting for the attraction of the mass 
directly above or below. 

The technique of using the integral of Green's functions 
rather than Green's functions themselves does not limit itself 
to gravity calculations only. The same technique can be used 
for all effects such as displacement, tilt, and strain calcula- 
tions to remove one singularity. It is not only limited to 
global tide models, but can be used regionally if the regional 
representations are given as areas of constant amplitude and 
phase. 

All ocean tidal contributions calculated in this study omitted 
the degree zero term (n=0) in order to impose mass conservation. 
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Because some of the observations were taken near the ocean, a 
small nonzero initial central angle was used. 

M2 OCEAN TIDE 

The global 1° square representation of the M2 constituent was 
obtained from E. W. Scbwiderski, of the Naval Surface Weapons 
Center, Dahlgren. This M2 model was generated by "hydrodynami- 
cal interpolation" (Schwiderski 1978). That is, in solving the 
Laplace tidal equations, more than 2,000 empirical tide gage ob- 
servations from continental and island stations were used to 
constrain the solution height amplitude and phase inside the 
grid compartments where tide gage observations were available. 
This feature is very important for studying ocean loading. Fre- 
quently, modification of global tide models is undertaken using 
more realistic local models because gravity measurements taken 
near coasts are sensitive to the local tide. However, for the 
results quoted here, no modifications of the Schwiderski M2 
model were made. The fine mesh size (l°xl°) was also an impor- 
tant consideration in choosing this particular model for this 
study. 


RESULTS 

Tidal Observations 

Normally, analyses of tidal gravity series are given in terms 
of amplitude factor 6 and phase tjj. The amplitude factor is the 
ratio of actual tidal response amplitude to the theoretical 
gravity value for a rigid Earth. Solutions of Earth modeling 
yield an amplitude factor of 6=1.16 and phase t|i=0 0 . The M2 
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signal in microgals at any time, t, can be generated from the 
expression 


fivl f)8 

Ag M2 (<M,t) = -6 v5/ .;. 6n g c M2 co S 2(«},) (9) 

• cos (n • p + 2k + i|i ) 

where p was defined earlier, K is east longitude, n=(2,0, 0,0, 0, 0) , 
* * , , 2 
g is magnitude of gravity in m/sec , a is the Earth semimajor 

axis in meters, and d> is latitude. The numerical value for 
is taken from Cartwright and Edden (1973). C M2 is equal to 
0.63192. Differing sign conventions for the local phase angle, 

»|), are found in the literature. For this reason (9) is expli- 
citly shown. When the phase angle, i|< , takes on negative values, 
it is regarded as a lag. 

The most precise tidal gravimeter results available are from 
the superconducting gravimeter studies of Warburton et al . 
(1975). These results are reported to have an accuracy of ±0.2% 
Tidal series with this gravimeter are available for La Jolla and 
Pinon Flat, Calif. 

Tidal gravimeter results for the Australian stations Alice 
Springs and Canberra were taken from Melchior (1978). Results 
for Bruxelles were taken from Melchior et al. (1976). Observa- 
tion results for Walferdange (Torge and Wenzel 1977), Potsdam 
(Altmann et al. 1977), and Mizusawa (Hosoyama 1977) were presen- 
ted at the Eighth International Symposium on Earth Tides, Bonn, 
September, 1977. 

The results of this study are given in table 1. The observa- 
tions in terms of amplitude factor and local phase are given in 
columns A. The observed values are reduced by the theoretical 
solid Earth values 6=1.16 and ijj=0 o and are given in columns B. 
Columns C show the predicted ocean tide contribution to the 
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gravity measurements using the Schwiderski M2 ocean tide model, 
the loading deformation coefficients of Farrell (1972) and 
Zschau (1978), and the integral Green's function technique pre- 
sented in this paper. By comparing columns B and C, one can see 
that the predicted ocean contribution to the gravity signal a- 
grees with the theoretical values (1.16,0°) to the 0.5-microgal 
level. These results are slightly better than the Australian 
results obtained by Bretreger and Mather (1978) using global 
ocean tide models of Hendershott and Zahel . 

Columns D represent corrected amplitude factors and phases 
under the assumption that the predicted ocean tide contributions 
do indeed properly model the ocean load. Except for Pinon Flat 
and Walferdange, the phases seem to show a negative trend. The 
magnitude of these phases is contrary to that predicted by 
Zschau (1978) for an imperfectly elastic Earth. His modeling 
shows that the lag in gravity measurements should be very small 
(order of 1/1,000°). These corrected results should not be 
taken too seriously, however. These phases represent measure- 
ment accuracies of 0.5 microgal or less, which is not the case. 
They may also be subject to common calibration errors. Further 
improvement is also possible in modeling the ocean tide in the 
open oceans where direct measurements of the ocean tidal ampli- 
tude and phase are sparse. 

Differences in Tidal Gravity 

Because of the quality of the observations at La Jolla and 
Pinon Flat, further investigation is indicated as a result of 
the disagreement between the observed tidal values at these lo- 
cations after correcting for the ocean contribution. The ocean 
predictions here seem to be slightly worse than those calculated 
by Warburton et al. (1975). Elimination of ocean tide effects 
from far afield is accomplished by subtracting the ocean effects 
between the two sets of observations. (La Jolla and Pinon Flat 
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are only 1° apart.) Table 2 shows the results of such differen- 
cing. One immediately notices that the differential predictions 
between Warburton et al. (1975) and the technique used here, 
along with the Schwiderski M2 model, are almost identical. 

Table 2. — M2 ocean tidal differences in gravity between 
La Jolla and Pinon Flat 



Amplitude 

Phase 


(pgai ) 

(degree) 

Observed 

2.25 

39.0 

Goad (this study) 

1.65 

38.3 

Warburton et al. (1975) 

1.61 

39.6 


Western European Gravity Comparisons 

Many tidal observations have been taken in western Europe 
during the past several years . Many of these observations were 
published by Melchior et al . (1976). Since then, Souriau (1979) 
published a set of corrections to the Melchior, Kuo, Ducarme ob- 
servation set for the effects of the ocean tidal loading. His 
procedure was to use the Green's functions of Farrell (1972) in 
conjunction with ocean tidal information obtained from digitized 
cotidal charts at 0.25° spacing for the neighboring seas. Co- 
tidal charts of several investigators were used to model the 
large water bodies. Table 3 gives these comparisons for western 
Europe. Notice that the predictions by Souriau and those 
computed using the Schwiderski M2 model with the load deforma- 
tion coefficients of Zschau are very similar. The major dif- 
ferences occur at Bordeaux and Cambridge where the effects of 
the oceans are rather large. Obviously the 0.25° resolution of 
nearby seas was an important ingredient in the tidal corrections 


13 



14 


at these two sites. Nevertheless, predictions of this quality 
using a global M2 model without any modifications are very good. 

CONCLUSIONS 

Integrals of Green's functions are more suitable for the 
special case of tidal loading where the tide is represented by 
areas or cells of constant amplitude and phase. Their use 
allows for the inclusion of the height of the instrument above 
sea level. This method also directly uses sequences of load de- 
formation coefficients which is advantageous for making compari- 
sons. The 1° square Schwiderski M2 ocean tide model predicts 
accurately ocean load perturbations of Earth tidal gravity obser- 
vations . Global calculations in this paper seem to be good to 
the 0.5-microgal level. Hopefully, with improvements in surface 
and space techniques, we are approaching a period when accurate 
checks between strain, tilt, gravity, and displacement observa- 
tions will be possible for the prediction of effects of ocean 
tide models . 
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